Comparaison du MIC de matrices de spectres Données human Serum prétraitées par SOAP par deux méthodes:
Lecture des matrices des spectres et initiatilisation des variables utiles
# Choix du répertoire où sont les matrices de données
pathdir="~/Dropbox/Partage_OMICS_COURSE/Data_stat2340/Human_Serum"
# Definition des matrices à lire dans ce répertoire
data_used=c("HumSerum_full_PEPS.csv","HumSerum_Manual_500P.csv")
spectral_matrices_names=c("HumSerum_full_PEPS","HumSerum_Manual_500P")
nmatrices=length(data_used)
# Lecture des matrices de données.
## On les stocke dans spectral_matrices_list
## On remplace les noms des colonnes par des noms sans le X
spectral_matrices_list=list()
for(i in 1:nmatrices)
{
spectral_matrices_list[[spectral_matrices_names[i]]]=as.matrix(read.table(file.path(pathdir,data_used[i]),header=T,sep=";",row.names=1,dec="."))
dimnames(spectral_matrices_list[[i]])[[2]]=as.numeric(substring(dimnames(spectral_matrices_list[[i]])[[2]],first=2))
cat("\n Lecture de la matrice : ",spectral_matrices_names[i])
}
##
## Lecture de la matrice : HumSerum_full_PEPS
## Lecture de la matrice : HumSerum_Manual_500P
# On prend la matrice des spectres communs dans les 2 bases (si ce n'est pas le cas au départ), on eleve les blancs des noms
dimnames(spectral_matrices_list[[1]])[[1]]=gsub(" ","",dimnames(spectral_matrices_list[[1]])[[1]]) # on enleve des blancs des noms des spectres.
namespectra=dimnames(spectral_matrices_list[[1]])[[1]]
countspectra=matrix(nrow=nmatrices+1,ncol=2)
dimnames(countspectra)=list(c(spectral_matrices_names,"Common Spectra"),c("Nb of spectra","Nb Buckets"))
countspectra[1,1]=length(namespectra)
countspectra[1,2]=dim(spectral_matrices_list[[1]])[2]
for(i in 2:nmatrices)
{
dimnames(spectral_matrices_list[[i]])[[1]]=gsub(" ","",dimnames(spectral_matrices_list[[i]])[[1]])
newnames=dimnames(spectral_matrices_list[[i]])[[1]]
compnames=is.element(namespectra,newnames)
countspectra[i,1]=length(newnames)
countspectra[i,2]=dim(spectral_matrices_list[[i]])[2]
namespectra=namespectra[compnames]
}
countspectra[nmatrices+1,1]=length(namespectra)
pander(countspectra)
| HumSerum_full_PEPS |
32 |
500 |
| HumSerum_Manual_500P |
31 |
500 |
| Common Spectra |
31 |
NA |
# On extrait les bons spectres et on les met dans le même ordre
for(i in 1:nmatrices)
{
spectral_matrices_list[[i]]=spectral_matrices_list[[i]][namespectra,]
}
# Creation des variables qui définissent les groupes
# Pour les données human serum c'est le 5ième charactère du nom du spectre
groupes=as.numeric(substring(dimnames(spectral_matrices_list[[1]])[[1]], 5, 5))
names_spectra=substring(dimnames(spectral_matrices_list[[1]])[[1]],1,5)
repetition=as.numeric(substring(dimnames(spectral_matrices_list[[1]])[[1]], 2, 2))
# 3. on dessine les spectres
for(i in 1:nmatrices)
{
print(spectral_matrices_names[i])
Draw(spectral_matrices_list[[i]],type.draw="signal",subtype="stacked",num.stacked=4)
}
## [1] "HumSerum_full_PEPS"








## [1] "HumSerum_Manual_500P"








Etude de la répétabilité spectrale
Between and Within group inertia
# Impression des résultats des calculs d'inertie
cat("Inerties non standardisées")
## Inerties non standardisées
pander(list(Inerties_non_standardisees=Inertiavalues,Inerties_en_PC=Inertiapc,Inerties_within_groups=InertiaPerGroup))
PCA
Eigenvalues
PCA Score Plots
draw.scores12 = vector("list")
draw.scores34 = vector("list")
for(i in 1:nmatrices)
{
draw.scores12[[spectral_matrices_names[i]]]=MBXUCL::DrawScores(PCA.res[[i]], drawNames=TRUE, type.obj = "PCA",createWindow=FALSE, main = paste0("PCA score plot for ",spectral_matrices_names[i]),pch = groupes, col = groupes,axes =c(1,2))
draw.scores34[[spectral_matrices_names[i]]]=MBXUCL::DrawScores(PCA.res[[i]], drawNames=TRUE, type.obj = "PCA",createWindow=FALSE, main = paste0("PCA score plot for ",spectral_matrices_names[i]),pch = groupes,col = groupes, axes =c(3,4))
}
## $HumSerum_full_PEPS

##
## $HumSerum_Manual_500P

## $HumSerum_full_PEPS

##
## $HumSerum_Manual_500P

PCA Loading Plots
draw.loadings12 = vector("list")
draw.loadings34 = vector("list")
for(i in 1:nmatrices)
{
draw.loadings12[[spectral_matrices_names[i]]] = MBXUCL::DrawLoadings(PCA.res[[i]], type.obj = "PCA",createWindow=FALSE, main = paste0("PCA loadings plot for", spectral_matrices_names[i]),axes = c(1:2), loadingstype="l", num.stacked = 2)[[1]]
draw.loadings34[[spectral_matrices_names[i]]] = MBXUCL::DrawLoadings(PCA.res[[i]], type.obj = "PCA",createWindow=FALSE, main = paste0("PCA loadings plot for", spectral_matrices_names[i]),axes = c(3:4), loadingstype="l", num.stacked = 2)[[1]]
}
draw.loadings12
## $HumSerum_full_PEPS

##
## $HumSerum_Manual_500P

## $HumSerum_full_PEPS

##
## $HumSerum_Manual_500P

Unsupervised clustering (MIC)
nClust = length(unique(groupes)) # nombre de clusters à rechercher
clustres=matrix(0,nrow=nmatrices,ncol=8)
dimnames(clustres)=list(spectral_matrices_names,c("DunnW","DunnKM","DBW","DBKM","RandW","RandKM","AdjRandW","AdjRandKM"))
for(i in 1:nmatrices)
{
print(spectral_matrices_names[i])
ClustMIC.res = MBXUCL::ClustMIC(Intensities = spectral_matrices_list[[i]], nClust = nClust, Trcl = groupes, Dendr = TRUE)
clustres[i,]=as.numeric(ClustMIC.res[1,1:8])
}
## [1] "HumSerum_full_PEPS"

## [1] "HumSerum_Manual_500P"

Table continues below
| HumSerum_full_PEPS |
0.8246 |
0.8246 |
0.6776 |
0.6776 |
1 |
1 |
| HumSerum_Manual_500P |
0.9332 |
0.9332 |
0.7824 |
0.7824 |
0.7054 |
0.7054 |
| HumSerum_full_PEPS |
1 |
1 |
| HumSerum_Manual_500P |
0.2324 |
0.2324 |
PLS-DA
nLVPLSDA = 4 # nombre du composantes du PLSDA
nrep=nlevels(as.factor(groupes)) # Nombre de réponses
PLSDA.res=list()
perf.plsda.RMSEP=matrix(nrow=nmatrices,ncol=nrep)
dimnames(perf.plsda.RMSEP)=list(spectral_matrices_names,paste0("Y",1:nLVPLSDA))
perf.plsda.R2=perf.plsda.RMSEP
perf.plsda.Q2=perf.plsda.RMSEP
for(i in 1:nmatrices)
{
cat("PLS-DA pour ",spectral_matrices_names[i])
PLSDA.res [[i]]= PLSDA(x = spectral_matrices_list[[i]], y = groupes, nLV = nLVPLSDA, drawRMSEP = TRUE)
perf.plsda.RMSEP[i,] = PLSDA.res[[i]]$RMSEP
perf.plsda.R2[i,]=PLSDA.res[[i]]$R2
perf.plsda.Q2[i,]=PLSDA.res[[i]]$Q2
}
## PLS-DA pour HumSerum_full_PEPS

## PLS-DA pour HumSerum_Manual_500P

pander(list(RMSEP=perf.plsda.RMSEP,R2=perf.plsda.R2,Q2=perf.plsda.Q2))
RMSEP:
| HumSerum_full_PEPS |
0.1145 |
0.0792 |
0.1061 |
0.1117 |
| HumSerum_Manual_500P |
0.3953 |
0.2478 |
0.426 |
0.2985 |
R2:
| HumSerum_full_PEPS |
0.9728 |
0.9873 |
0.971 |
0.9662 |
| HumSerum_Manual_500P |
0.619 |
0.8414 |
0.6475 |
0.6968 |
Q2:
| HumSerum_full_PEPS |
0.9316 |
0.9672 |
0.9356 |
0.9349 |
| HumSerum_Manual_500P |
0.1839 |
0.6794 |
-0.0381 |
0.5347 |
Score plots de la PLS_DA
PLSDA.scores= vector("list")
for(i in 1:nmatrices)
{
PLSDA.scores[[spectral_matrices_names[i]]]=DrawScores(PLSDA.res[[i]], drawNames=TRUE, type.obj = c("PLSDA"),
createWindow=FALSE, main = paste0("PLSDA score plot for ", spectral_matrices_names[i]),pch = groupes,col = groupes, axes =c(1,2))
}
PLSDA.scores
## $HumSerum_full_PEPS

##
## $HumSerum_Manual_500P

Loading plots de la PLS_DA
PLSDA.loadings= vector("list")
for(i in 1:nmatrices)
{
PLSDA.loadings[[spectral_matrices_names[i]]]=DrawLoadings(PLSDA.res[[i]], type.obj = "PLSDA",
createWindow=FALSE, main = paste0("PLSDA loadings plot for", spectral_matrices_names[i]),
axes = c(1:2), loadingstype="l", num.stacked = 2)
}
PLSDA.loadings
## $HumSerum_full_PEPS
## $HumSerum_full_PEPS[[1]]

##
##
## $HumSerum_Manual_500P
## $HumSerum_Manual_500P[[1]]
